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Abstract:  A  model  of  transient  modal  instability  in  fiber  amplifiers  is 
presented.  This  model  combines  an  optical  beam  propagation  method  that 
incorporates  laser  gain  through  local  solution  of  the  rate  equations  and 
refractive  index  perturbations  caused  by  the  thermo-optic  effect  with  a  time- 
dependent  thermal  solver  with  a  quantum  defect  heating  source  term.  This 
model  predicts  modal  instability  in  a  285  Watt  fiber  amplifier  characterized 
by  power  coupling  to  un-seeded  modes,  the  presence  of  stable  and  unstable 
regions  within  the  fiber,  and  rapid  intensity  variations  along  the  fiber. 
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1,  introduction 

Recent  evidence  suggests  that  dynamic  thermal  modal  instability  (TMI)  affects  the  beam 
quality  of  large  mode  area  high  average  power  fiber  amplifiers  [1,23,4,5].  The  resulting 
degradation  severely  limits  the  usefulness  of  these  devices.  Efforts  to  improve  fiber  amplifier 
designs  to  mitigate  this  degradation  benefit  from  numerical  models  that  reveal  the  effects  of 
amplifier  parameters  on  performance.  Several  previous  investigations  have  all  but  isolated  the 
cause  to  refractive  index  gratings  caused  by  spatial -temporal  temperature  variations  that  arise 
from  quantum  defect  heating  through  the  thermo-optic  effect  [6,7*8,9,10]*  A  key  aspect  of  the 
physical  description  of  this  phenomenon  is  the  time  behavior  of  the  temperature  profile  within 
the  fiber*  Under  the  assumption  that  the  temperature  at  any  point  within  the  fiber  is  a  periodic 
function  of  time,  a  steady-state  solution  to  the  optical  wave  equation  arises  that  describes  the 
process  of  stimulated  thermal  Rayleigh  scattering  (STRS)*  Numerical  models  have  been 
presented  that  capture  the  essential  features  of  STRS  in  fiber  amplifiers  [9,10,11,12], 
Observations  of  periodic  modal  fluctuations  in  the  output  beams  of  amplifiers  operating  near 
the  TMI  threshold  are  consistent  w  ith  the  STRS  picture  of  TMI  [3,7].  This  regime  of  TMI 
rehes  on  the  amplifier  optical  and  temperature  fields  achieving  a  dynamic  equilibrium  which 
may  be  unstable  with  respect  to  external  perturbation  [9]. 

Although  significant  understanding  is  provided  by  the  steady-state  model,  several  factors 
motivate  further  study  of  the  transient  behavior  of  TMI.  Amplifiers  operated  with  a  rapid 
turn-on  may  not  reach  steady-state  equilibrium  prior  to  the  onset  of  TMI  thus  the  threshold 
calculated  using  the  periodic  method  may  be  higher  than  values  that  characterize  realistic 
operating  conditions.  One  definition  of  the  TMI  threshold  has  been  given  in  terms  of  the 
buildup  time  required  to  observe  the  effect  [5].  Also,  observed  TMI  is  often  characterized  by 
chaotic  behavior  [3,4,7]  not  able  to  be  described  in  the  context  of  a  periodic  model.  The  goal 
of  the  work  presented  here  is  to  present  a  computational  model  that  can  enable  a  greater 
understanding  of  the  chaotic  and  transient  regimes  of  TMI*  This  model  incorporates  a  3D 
beam  propagation  method  {8 PM),  including  laser  gain,  for  the  spatial  evolution  of  the  signal 
field  as  well  as  a  3D  time-dependent  thermal  model.  This  model  is  also  capable  of  describing 
the  steady-state  periodic  case. 

This  paper  begins  with  a  description  of  the  BPM  which  is  of  a  new  type,  followed  by  a 
description  of  the  thermal  solver,  also  of  a  new  type*  Example  results  are  then  presented  and 
discussed* 

2*  Hybrid  finite  element- harmonic  beam  propagation  method 

The  problem  at  hand  requires  the  beam  propagation  calculation  to  be  performed  over  the 
entire  length  of  the  fiber  thousands  of  times,  once  for  each  time  step.  Each  propagation  step 
must  occur  over  a  fraction  of  the  inter- mo  dal  beat  length  of  the  fiber  leading  to  a  requirement 
to  evaluate  upwards  of  one  million  BPM  steps*  Thus  the  computational  speed  of  the  BPM  is 


at  a  premium.  Recently,  a  new  type  of  BPM  was  presented,  the  azimuthal  harmonic 
expansion  beam  propagation  method  [13],  which  enables  much  faster  computation  than  either 
the  split-step  fast  Fourier  transform  BPM  (FFT-BPM)  or  finite  difference  BPM  (FD-BPM). 
Ihis  method  is  based  on  a  hybrid  discretization  in  cylindrical  coordinates  that  combines  a 
Finite-difference  discretization  in  the  radial  coordinate  with  a  harmonic  Fourier  expansion  in 
the  azimuthal  angle.  However,  a  ID  finite  element  approach  in  the  radial  direction  aids 
treatment  of  fiber  designs  with  discontinuous  refractive  index  profiles,  such  as  photonic 
crystal  fibers,  as  well  as  enabling  variable  grid  spacing  to  conserve  the  total  number  of 
unknowns  in  a  natural  fashion.  A  simplification  of  such  a  method  has  been  previously  applied 
to  the  study  of  stimulated  Brillouin  scattering  in  fibers  [14]. 

The  starting  point  is  the  scalar  optical  wave  equation  Eq.  (1)  valid  for  weakly-guiding 
fibers 

V-E(r)+n2  (r)tfE(r)  =  0,  (1) 

where  k(t  -wfc  is  the  free  space  wave-vector  and  n(r)i$  the  refractive  index  distribution. 
Invoking  the  slowly-varying  envelope  and  paraxial  approximations  yields  Eq.  (2) 

-2'7?J;  +  V,3  +n2(r,<p,z)k2 -02 

where  r)  now  refers  to  the  slowly-varying  envelope,  ft  is  the  wave -vector  describing 

rapid  oscillations  in  the  propagation  direction,  and  V^is  the  transverse  Laplaeian  operator. 

Setting  dE/  cz  =  %dnidz  -  0  yields  an  equation  Eq,  (3)  for  the  stationary  propagating  modes 
of  the  fiber 


E{r,<p,z)=  0  (2) 


[V;  +/r  {>\<p)kl  ]£(>-,<?)  =  0.  (3) 

To  obtain  a  finite  element  solution  of  this  equation  the  variational  form  Eq.  (4)  is  introduced 


5  =  ^JJ[“I V<£(r’^r  +”20'*^)Ao  \E(r><pf  -01 

where  Q  represents  the  fiber  cross-section.  Expanding  as  shown  in  Eq.  (5) 

Q 

E(r,<p)  =  X  £,(r)exP['?d 
where  O  is  the  truncation  order  of  the  series  leads  to  Eq.  (6) 
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where  Eq.  (7)  and  Eq.  (8)  define  the  refractive  index  in  Eq.  (4), 
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wrhere  is  a  suitable  background  refractive  index.  With  these  definitions  the  variational 
form  is  given  by  Eq.  (9) 


(9) 


To  accomplish  the  radial  integral  it  is  convenient  to  divide  the  integration  domain  into  ID 
finite  elements  represented  by  N  -1  line  segments  with  end  points  rk  :  k  <s  1  yN  .  The  integrals 
over  the  fiber  cross  section  are  then  approximated  as  shown  in  Eq.  (10)  and  Eq.  (11) 

JJ/(>,)exp[Mt^  ~  +2~--j(^i  ~rt )  o°) 

Him- exp \iq<p]dA  »  7rS(q)  jji  I  (rw  -r*  )  (11) 

ft  k=\  \  '*  +  1  fk  ) 

The  variational  form  is  then  conveniently  expressed  using  matrix  notation  as  shown  in  Eq. 

(12) 


S  =  -  ~E*  -[k  -/Tm]-  E  (12) 

the  minimi  zation  of  which  leads  to  the  generalized  eigenproblem  given  in  Eq,  (13) 

Fk-^3!M]-E  =  0  (13) 

where  Eis  a  vector  of  dimension  (2£)+1)jV  and  Kand  M  are  square  matrices  with  the  same 

dimension.  If  now  the  envelope  is  allowed  to  vary  slowly  along  the  length  of  the  fiber  such 
that  BE/dz^  0,  then  the  longitudinal  wave  vector  value  k  is  freely  chosen  leading  to  the 
matrix  form  of  the  paraxial  wave  equation  as  shown  in  Eq.  (14). 

-2/£M  —  +  K  -  A2M  j  •  E  =  0  (14) 

Differencing  E with  respect  to  Zand  taking  one  explicit  and  one  implicit  step  each  of 
distance  Az/2  results  in  the  usual  Crank-Nicholson  update  rule  for  propagation  given  in  Eq. 
(15) 

(M-^(Km-A2M))Em  (15) 

where  the  cumulative  propagation  distance  is  z  =  (m-  l)Ar  ;  m  e  l,M  .  In  this  scheme,  the 

waveguide  refractive  index  profile,  the  gain  within  the  core,  and  the  thermally- induced 
refractive  index  profile  are  incorporated  as  shown  in  Eq,  (16). 

8rf  (/%  <p,z)  =  (r,<p)  +  inug  ( r ,  <p,z)  +  2n0  ~  A  T(r,  <p,z )  (16) 

where  is  the  waveguiding  index  profile,  g  is  the  optical  gain,  drt/dT  is  the  thermo¬ 
optic  coefficient,  and  AT  is  the  local  temperature  rise.  To  a  good  approximation,  the  gain  is 
instantaneously  determined  by  the  local  pump  and  signal  intensities  while  the  local 
temperature  rise  depends  on  the  past  history  of  heat  deposited  in  the  neighborhood  of  each 
point  within  the  fiber  according  to  the  time-dependent  heat  equation. 

The  matrix  K  changes  at  each  step  due  to  its  dependence  on  the  signal  intensity  profile. 
To  keep  the  gain  contribution  space-centered  during  the  propagation  step,  the  field  is  updated 
first  using  the  initial  matrix  Km  for  both  the  explicit  and  implicit  steps.  Then  these  field 
values  are  used  to  calculate  Km+i  which  is  then  used  to  repeat  the  implicit  step.  Furthermore, 
the  average  refractive  index  change  can  increase  significantly  over  the  length  of  the  fiber 


necessitating  resetting  the  reference  refractive  index  between  steps.  This  is  done 
automatically  in  the  computer  implementation. 

3*  Gain  and  heat  deposition 

The  following  discussion  pertains  to  active  media  for  which  the  atomic  energy  level 
populations  are  determined  by  two-level  rate  equations  where  instantaneous  decay  from  the 
lower  level  to  the  ground  state  is  assumed.  Under  these  conditions  the  local  fraction  of  ions  in 
the  upper  state  is  given  by  Eq,  (17) 


/  er 


hv. 


£  +  W 


hv„ 


I  ^  +  ^,v ) T  ^ p  (*7^  +  ^ ) T 

hVr  hv„ 


(17) 


where  Isp  are  the  signal  and  pump  intensities,  hvsp  are  the  energies  of  each  signal  and  pump 

photon,  a  are  the  absorption  and  emission  cross-sections  at  the  pump  and  signal  wavelengths, 
and  t  is  the  upper-state  lifetimes.  The  gain  is  then  given  by  Eq.  (18) 

g(t\<p,z)  =  Nim \f{r,<p,z )cr,s  -(l-/(r,^r)cru)]  (18) 


where  NltXi  is  the  doping  density.  Similarly,  the  pump  gain  is  given  by  Eq.  (19) 


gp  (r,<p,z)  =  Nim  [/(< r ,v,z)trv  -(l  -/(r,^z)o-„,)].  (19) 

For  cladding-pumped  amplifiers,  the  paraxial  approximation  is  generally  not  valid  for  the 
pump  light  guided  in  the  cladding  due  to  its  large  numerical  aperture.  Therefore,  constant 
pump  intensity  throughout  the  cladding  is  assumed  to  be  maintained  through  the  process  of 
mode  mixing  so  that  the  evolution  of  the  pump  intensity  is  given  by  Eq,  (20) 

47/.(z)  =  ±X -HSr(r'V'x)Ip(r'VtZ)di  (20) 

where  Ac]ad  is  the  area  of  the  pump  cladding  and  the  ±  represents  counter  versus  co-pumping 
arrangements.  The  heat  load  due  to  the  quantum  defect  is  then  given  by  the  number  of  ions 
going  through  the  cycle  multiplied  by  the  rate  that  they  go  through  the  cycle  and  the  quantum 
defect  for  each  photon  as  shown  in  Eq,  (21) 

^--lj  f{r,<p,z)<rjs(r,<p,z)  (21) 


Q{r  ,<p,z)  =  N 


The  normalization  of  the  optical  fields  is  chosen  so  that  L  =  |£|2 .  Starting  from  some  known 

temperature  distribution,  the  temperature  evolves  in  time  at  each  spatial  point  according  to  the 
time -dependent  heat  equation  given  by  Eq.  (22) 

^T(r,p,z,t)  =  -^-V2T(r,<p,z,t)+~Q(r,p,z,f)  (22) 

vi  pL  pL 

where  K  is  the  thermal  conductivity,  p  is  the  mass  density,  and  C  is  the  heat  capacity  within 
the  fiber.  This  equation  is  amenable  to  a  solution  approach  similar  to  that  used  for  the  optical 
propagation  equation.  Both  T  and  O  are  approximated  in  rand  <p  using  the  ID  finite  element 
scheme  in  the  radial  direction  and  the  azimuthal  harmonic  expansion.  The  longitudinal 
component  of  the  Laplacian  is  evaluated  using  finite  differences  in  the  temperature  values  at 
the  z  points  defining  the  optical  beam  propagation  steps  as  shown  in  Eq.  (23) 
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To  a  good  approximation,  the  thermal  properties  of  the  fiber  are  constant  throughout  the 
cross-section.  Significantly,  this  causes  the  different  azimuthal  harmonic  terms  to  de-couple 
yielding  (20  +  1)  independent  Crank-Nicholson  time  step  update  equations  tor  the 
temperature  distribution  as  shown  in  Eq.  {24} 


At  K 

Mth  +  2  pf  - 
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(24) 

where  l  Is  the  time  step  index.  At  is  the  time  step,  q  is  the  harmonic  component,  T  ,  are 
vectors  of  length  A/  A  and  MT1|,K?THare  sparse  matrices  of  this  same  dimension.  These 
sparse  matrices  are  assembled  by  inserting  blocks  of  values  that  when  multiplied  by  the 
approximate  the  integrals  over  the  fiber  cross-section  as  described  by  Eq.  (25)  and  Eq.  (26): 
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where  AAk  is  the  area  of  the  annulus  with  inner  and  outer  radii  of  rk  and  rt+l ,  vkm  are  column 

vectors  that  have  two  non-zero  entries  each  with  value  I  /  2  at  locations  in  *  N  +  k  and 
m*N+k  + 1 ,  and  uA.  m  are  column  vectors  that  have  two  non -zero  entries  with  values  -1  and  I 

at  locations  mN  +  km\d  m- N +k  +  \ .  The  choice  to  order  the  vectors  X, ,  first  in  k  is 

arbitrary  and  leads  to  contiguous  blocks  non-zero  entries  of  dimension  N  in  the  matrices. 

To  keep  the  time  step  update  centered  in  time,  the  harmonic  components  of  the  heat  load 
given  by  Eq,  (27) 
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must  be  evaluated  at  the  beginning  and  end  of  each  time  step.  Therefore  the  updated 
temperature  is  calculated  first  using  the  initial  heat  distribution.  This  temperature  distribution 
is  then  used  to  update  the  optical  field  leading  to  the  updated  heat  distribution.  This  is  then 
used  to  repeat  the  temperature  update,  this  time  incorporating  the  heat  distributions  at  the 
beginning  and  end  of  the  time  step.  To  start  the  simulation  an  initial  temperature  distribution, 
input  mode  field,  and  input  pump  power  must  be  specified. 

The  input  mode  field  is  taken  to  be  a  superposition  of  guided  modes  with  pre- specified 
amplitudes.  Additionally,  a  time-dependent  phase  shift  may  be  applied  to  one  or  more  of  the 
input  modes  to  incorporate  a  frequency  shift  between  modes  as  shown  in  Eq.  (28) 

E(r,<p,0,t)  =  ZV^£»(r’<p)exp| jf»  (0] 


(28) 


where  Pn  and  Et}  (r,  <p)  are  the  power  and  mode  field  distribution  of  guided  mode  n  ,  This  has 

been  shown  to  induce  inter-modal  coupling  through  the  process  of  stimulated  thermal 
Rayleigh  scattering  [9*10*1 1,12]. 

4.  Computational  implementation 

The  model  described  above  is  a  3+ ID  fiber  amplifier  model  with  a  hybrid  discretization 
scheme  in  cylindrical  coordinates  employing  finite  elements  in  the  radial  direction,  a 
harmonic  expansion  in  the  azimuthal  direction,  and  a  finite  difference  grid  in  the  direction  of 
optical  propagation.  For  a  fiber  length  of  1.63  meters  with  a  beat  length  of  223  millimeters, 
100  samples  per  beat  length  lead  to  a  longitudinal  grid  with  7300  points.  The  finite  element 
approach  in  the  radial  dimension  permits  variable  radial  point  spacing  thus  reducing  the 
number  of  radial  points  to  about  100  that  span  the  entire  fiber  cladding  radius,  resulting  in  a 
matrix  order  of  approximately  7,3  *105  for  each  of  the  (2£?+l)  thermal  matrices.  Since  the 

region  involved  in  optical  propagation  is  only  at  the  center  of  the  fiber,  a  subset  of  the  thermal 
radial  points  is  used  in  the  optical  propagation  equations.  This  region  is  typically  2-3  times 
the  core  size.  Fig,  1  depicts  the  computational  steps  required  to  implement  this  method. 


Fig.  I .  Sequence  oF  required  to  carry  out  the  time  dependent  amplifier  simulation. 

This  approach  can  be  used  for  any  waveguide  profile  for  which  the  harmonic  expansion 
can  be  calculated,  however,  those  with  stronger  azimuthal  variations  require  higher  truncation 
orders.  It  is  convenient  to  use  separate  truncation  orders  Q  for  the  optical  and  thermal 
problems.  For  modeling  photonic  crystal  fibers  with  small  capillaries  truncation  orders  of 
0  □  10-20  are  appropriate  while  for  step  index  fibers,  Qop]  can  be  chosen  to  include  the 
highest  azimuthal  order  propagating  modes  supported  by  the  fiber.  For  typical  large  mode 
area  fibers  QtJV3  0  2“4  suffices,  Likewise,  the  averaging  effect  resulting  from  heat  diffusion 


means  that  Qlhum  □  2-4  for  all  types  of  fibers.  Nevertheless,  the  orders  may  be  increased 
until  the  desired  level  of  convergence  is  achieved.  The  simulation  time  is  most  sensitive  to 
£?opt  because  this  affects  the  speed  of  the  beam  propagation  portion  of  the  calculation  which  is 

the  most  time  consuming.  The  simulation  time  of  a  large  pitch  photonic  crystal  fiber  was 
observed  to  be  about  five  times  greater  than  that  of  a  step  index  fiber. 

Within  this  approach,  the  primary  computational  tasks  are  sparse  matrix  multi  plication  and 
sparse  linear  system  solution.  Various  software  packages  are  available  for  these  tasks  that 
take  advantage  of  modern  multi -processor  high  performance  computing  architectures.  The 
implementation  reported  here  relies  on  parallelization  using  the  message  passing  interface 
(MP1)  [15]*  Within  MPI,  each  collective  task  on  the  processor  grid  is  mediated  through  an 
MP1  communicator.  To  solve  all  harmonic  components  of  the  thermal  problem 
simultaneously,  the  processor  grid  is  divided  into  sub-communicators,  one  for  each  harmonic 
component.  Therefore,  as  long  as  processors  are  available,  the  thermal  harmonic  order  may 
be  increased  with  negligible  increase  in  the  time  required  to  accomplish  the  thermal  update. 

During  the  beam  propagation  portion  of  the  calculation,  the  temperature  distribution 
remains  stationary  in  memory  scattered  across  the  processing  grid  and  the  required 
temperature  values  are  broadcast  to  all  processes  at  each  spatial  propagation  step.  The  hybrid 
discretization  scheme  greatly  reduces  the  number  of  floating  point  operations  required  for 
each  update  step  compared  to  other  discretization  techniques.  Therefore  the  optical  solution  is 
obtained  faster  by  the  group  of  processors  on  each  sub-communicator  performing  its  own 
optical  update  than  by  spreading  the  optical  update  matrix  operations  across  the  entire 
processing  grid  due  to  the  communications  overhead  required.  As  the  thermal  order  is 
increased,  additional  communications  time  is  required  to  retrieve  the  additional  temperature 
information  for  the  optical  update  step  leading  to  a  gradual  increase  in  the  total  time  to 
solution, 

5*  Example  results 

As  a  first  example,  it  makes  sense  to  compare  the  results  of  this  model  to  a  prior  mode! 
that  employed  coupled  mode  theory  for  the  optical  propagation  [7],  The  double-clad 
Ytterbium-doped  amplifier  parameters  are  given  in  Table  1.  Due  to  the  calculation  speedup,  a 
step  index  fiber  approximately  equivalent  to  the  photonic  crystal  fiber  discussed  before  was 
simulated. 


Table  1,  Simulated  Fiber  Parameters 


Parameter 

Value 

core  diameter 

74  jim 

pump  cladding  diameter 

170  Jim 

outer  cladding  diameter 

400  Jim 

core  numerical  aperture 

0.03 

fundamental  mode  field  area 

2750  jim3 

beat  length  LPoi-^Pn 

22 3  mm 

fiber  length 

1.63  m 

yb+5  doping  concentration 

3.5 xlO25  m'3 

signal  wavelength 

1 .064  pm 

pump  wavelength 

0.977  pm 

signal  power  LPot 

9,5  W 

signal  power  LPn 

0,5  W 

pump  power 

357  W 

signal  emission  cross-section 

3, 58*1  O'25  m3 

signal  absorption  cross-section 

6.00x1  O'37  m2 

pump  emission  cross-section 

1.87X10'24  m3 

pump  absorption  cross -sect ion 

1.53x10'34  m3 

upper  state  lifetime 
thermal  conductivity 
heat  capacity 
mass  density 
thermo -optic  coefficient 
heat  sink  temperature 
z  grid  spacing 
time  step 

total  simulation  time 
thermal  radial  points 
optical  radial  points 
doped  core  radial  points 
optical  azimuthal  order 
thermal  azimuthal  order 


850  ps 
1.38  W/m-K 
703  J/kg-K 
2200  kg/m3 
].2xl(r 
300  K 
2+23x10"4m 
L0*10'5s 
10  ms 
101 
56 
30 
3 
2 


Furthermore,  there  was  no  material,  scattering,  or  bend  loss  assumed.  The  outer  boundary  of 
the  fiber  cladding  is  maintained  at  a  fixed  reference  temperature  thus  assuming  perfect 
conductive  cooling.  Also,  the  pump  and  signal  line  widths  were  assumed  negligible. 

A  counter-pumped  configuration  was  assumed  for  this  simulation  in  which  the  total  pump 
absorption  throughout  the  fiber  was  approximated  to  be  constant  so  that  the  pump  intensity  at 
the  seeded  end  could  be  set  as  a  boundary  condition. 
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Fig.  2,  Plot  of  the  optical  intensity  as  a  function  of  length  at  the  center  of  the  fiber  core  and  at 
points  offset  from  the  center  of  the  fiber  core  by  half  tire  core  radius  after  a  simulation  lime  of 
10  ms.  The  accompanying  movie  (Media  1 )  shows  the  evolution  of  these  intensity  probes  over 
a  time  period  of  10  ms  The  inset  shows  the  optical  field  at  the  output. 

Fig.  2  and  the  accompanying  movie  show  the  time  dependence  of  the  optical  intensity  within 
the  fiber  core  and  the  output  intensity  distribution  over  time.  The  simulation  was  started  with 
all  of  the  power  launched  into  the  LPDI  mode  and  then  the  LPn  was  gradually  introduced  over 
a  period  of  0. 1  ms.  The  introduction  of  the  higher  order  mode  immediately  causes  the  output 
to  become  unstable.  In  agreement  with  the  coupled  mode  theory  results  [7]  the  instability  is 
confined  to  the  latter  portion  of  the  amplifier  and  is  more  severe  near  the  output  end. 
Furthermore,  as  time  progresses  the  stable  region  grows  in  extent  as  the  thermally  induced 
grating  reaches  stable  equilibrium.  This  equilibrium  region  is  possible  due  to  the  fact  that  the 
two  modes  are  launched  at  the  same  optical  frequency  and  so  no  moving  grating  exists  at  the 
seed  end  [6],  This  is  also  in  agreement  with  the  coupled  mode  theory  results  [7]. 

While  there  are  some  similarities  between  the  appearances  of  Fig.  I  and  Fig.  2  of  [7] 
different  quantities  are  being  plotted.  Fig.  I  shows  the  local  intensity  probe  at  three  points 


while  the  prior  published  figure  shows  the  modal  content.  The  full  optical  field  that  is 
required  to  calculate  the  modal  content  was  not  stored  at  every  longitudinal  position  for  every 
time  step  due  to  the  large  file  size  that  would  be  required.  Nevertheless*  rapid  spatial 
oscillations  characteristic  of  non-adiabatic  power  transfer  [8]  are  evident  in  both.  The  beam 
propagation  model  described  here  also  captures  the  effect  of  thermal  lensing  on  the  mode  field 
area  of  the  fundamental  mode  of  the  fiber.  The  mode  field  area  of  the  cold  fiber  is  2750  pm3 
and  at  the  output  end  of  the  heated  fiber,  the  area  decreased  to  1960  pm2. 
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Fig.  3.  Modal  content  of  the  amplifier  output  over  time  as  well  as  total  power  from  the 
integrated  optical  output  intensity  and  the  sum  of  the  first  8  modal  powers. 

Fig.  3  shows  the  time  evolution  of  the  output  modal  decomposition,  it  is  significant  to 
note  that  even  though  only  two  transverse  modes  are  launched  at  the  input,  power  scatters  into 
other  modes  almost  immediately.  Since  the  thermal  boundary  condition  was  symmetric  in 
this  simulation,  one  would  expect  that  the  output  would  retain  mirror  symmetry  about  the  axis 
of  the  launched  LPn  mode  but  this  was  not  the  case.  This  was  attributed  to  the  fact  that  the 
sparse  linear  solvers  used  to  accomplish  the  beam  propagation  were  of  the  iterative  variety  so 
that  each  propagation  step  and  thermal  update  introduces  error  into  the  optical  fields  and 
temperature  profile  depending  on  the  chosen  convergence  tolerance.  It  is  tempting  to 
immediately  propose  a  correspondence  between  residual  error  in  the  equation  solutions  and 
unavoidable  imperfections  in  the  optical  waveguide  structure.  A  fundamental  question 
regarding  the  observed  instability  is  whether  the  symmetric  solution  is  unstable  with  respect  to 
power  present  in  non -symmetric  modes  and  the  results  presented  here  suggest  that  it  is.  This 
matter  warrants  further  study  that  is  beyond  the  scope  of  this  paper. 

The  total  output  power  was  found  to  fluctuate  in  disagreement  with  some  reports  [4].  This 
is  attributable  to  the  fact  that  the  spatial  fluctuations  of  intensity  within  the  fiber  inhibit 
efficient  extraction  of  the  available  gain  thus  causing  the  uniform  pump  absorption 
approximation  to  break  down.  This  observation  suggests  an  additional  mechanism  for 
feedback  that  can  drive  the  instability.  It  has  been  shown  that  amplitude  fluctuations  at  the 
seeded  end  of  amplifiers  can  cause  a  steady -state  amplifier  solution  to  become  unstable  [9].  If 
fluctuations  in  the  spatial  intensity  cause  varying  amounts  of  pump  to  be  absorbed  in  a 
counter-pumped  configuration,  this  would  effectively  modulate  the  seed  level  present  near  the 
input  end  of  the  amplifier  thus  providing  the  feedback  mechanism  modulating  the  amplitude 
causing  instability. 

It  is  also  apparent  that  the  incorporation  of  the  beam  propagation  model  significantly 
lowers  the  instability  threshold.  The  pump  powrer  of  357  W  in  the  simulation  here  was  below 
the  previously  observed  approximate  instability  threshold  of  1060  W  [7]  for  the  conduct! vely- 
cooied  case  and  yet  instability  is  clearly  present.  This  suggests  that  models  of  the  transient 


regime  based  on  coupled  mode  theory  should  include  all  guided  modes  in  order  to  accurately 
capture  instability  behavior.  This  also  agrees  with  the  observation  that  larger  cores  that 
support  more  transverse  modes  exhibit  lower  instability  threshold  powers.  Experimentally- 
observed  instability  thresholds  for  this  type  of  amplifier  are  yet  lower  so  clearly  the 
incorporation  of  the  beam  propagation  model  improves  prospects  for  agreement  of  the  model 
with  observed  thresholds. 


Fig.  4.  Frequency  spectra  of  the  first  five  modes  of  the  simulated  fiber. 

Fig.  4.  shows  the  frequency  spectra  of  the  transverse  modes  at  the  fiber  output.  No  sharp 
frequency  peaks  are  visible  indicating  a  degree  of  randomness  in  the  output.  Also,  the  high 
frequency  tail  was  more  prominent  in  the  lower-order  modes,  which  is  not  immediately 
evident  just  by  examining  the  time  series  shown  in  Fig.  3. 


Fig.  S.  Temperature  distribution  as  a  function  of  length  at  the  fiber  core  and  at  points  offset 
from  the  center  of  the  fiber  core  by  half  the  core  radius.  The  accompanying  movie  (Media  2) 
shows  the  evolution  of  these  temperature  probes  over  a  time  period  of  10  ms. 

The  temperature  was  recorded  in  lime  at  the  same  sampling  points  as  the  optical  intensity 
throughout  the  fiber  and  is  presented  in  Fig,  5,  As  expected,  the  temperature  is  static  in  the 
region  of  the  fiber  not  experiencing  instability  but  is  irregular  where  the  optical  intensity 
fluctuates.  Due  to  the  length  of  the  thermal  diffusion  time,  the  temperature  cannot  keep  up 


with  the  rapidly  fluctuating,  intensity.  Rather,  it  exhibits  slowly  varying  behavior  that  captures 
a  windowed  time  average  of  the  quantum  defect  heating  profile.  Both  the  center  temperature 
and  the  amplitude  of  the  temperature  fluctuations  are  irregular  and  not  uniform  across  the  core 
as  revealed  by  the  difference  between  the  profiles  at  the  core  center  and  off  center.  Observing 
the  time  evolution  of  the  temperature  profile  from  the  beginning  of  the  simulation  where  all  of 
the  power  was  in  the  LP^h  mode  reveals  an  immediate  overall  temperature  drop  as  power  is 
first  introduced  into  the  LPn  mode  and  then  coupled  into  the  other  higher  order  modes.  This 
is  consistent  with  the  prior  observation  that  higher  order  modes  create  a  lower  heat  load  than 
the  fundamental  [7]  and  therefore  transfer  of  power  out  of  the  fundamental  mode  decreases 
the  optical  length  of  the  amplifier. 

6*  Model  limitations,  and  future  investigations 

The  model  presented  here  has  some  limitations  in  its  current  form  that  deserve  discussion. 
The  first  is  that  bending  losses  and  bend-induced  mode  distortion  are  not  accounted  for.  An 
effective  index  gradient  may  be  incorporated  in  the  usual  fashion  to  account  for  bending  [16]. 
This  breaks  the  azimuthal  symmetry  of  the  step-index  profile  leading  to  additional  terms  in 
the  harmonic  azimuthal  expansion  of  the  index.  The  optical  boundary  condition  use  here  is 
lossless  and  perfectly  reflecting.  It  is  possible  that  transparent  boundary  conditions  [17]  could 
be  adapted  to  the  paraxial  scalar  beam  propagation  presented  here  as  was  done  for  a  previous 
similar  method  [13].  Furthermore,  material  loss  is  not  accounted  tor.  In  the  current 
implementation,  the  doped  region  of  the  core  must  be  assumed  to  be  uniform  and  circular. 
While  this  is  true  of  some  step  index  large  mode  area  fibers,  micro-structured  fibers  often 
have  some  additional  structure  that  arises  from  the  stack  and  draw  fiber  fabrication  method. 
As  was  already  mentioned  above,  the  presented  model  does  not  hold  the  launched  pump 
perfectly  constant  for  counter-pumped  configurations  which  are  practically  a  necessity  for 
minimizing  the  effective  amplifier  length  to  maximize  non-linear  thresholds.  An  efficient 
method  for  overcoming  this  limitation  in  the  transient  regime  is  not  immediately  apparent 
Finally,  in  its  current  form,  this  model  requires  high-performance  computing  resources  to  use. 
It  is  not  suitable  for  desktop  computers. 

The  model  presented  here  can  enable  numerous  future  investigations  a  few  of  which  are 
briefly  discussed  below.  First  it  is  important  to  verify  that  in  the  absence  of  noise,  quantum 
or  otherwise,  the  amplifier  can  reach  stable  thermo-optic  equilibrium.  While  a  previous 
investigation  predicted  dynamic  instability  even  in  the  absence  of  such  noise  [7],  this  could 
have  been  due  to  limits  in  the  achievable  accuracy  of  the  solutions  of  the  derived  equations. 
This  instability  was  also  predicted  to  occur  at  power  levels  much  higher  than  those  at  which 
has  been  observed  suggesting  that  noise  plays  an  important  role  in  the  origin  of  instability. 

This  model  could  also  be  used  to  analyze  the  effects  of  different  types  of  noise  on  the 
onset  of  instability.  This  could  include  frequency  offsets  between  modes,  fluctuations  of 
pump  and  seed  powers,  and  other  time  varying  launch  conditions.  Build  up  and  decay  times 
of  instability  could  also  be  studied.  The  potential  for  increasing  the  instability  threshold 
through  advanced  fiber  designs  could  also  be  studied.  For  example  the  model  can  treat  large 
pitch  photonic  crystal  fibers  that  rely  on  de- localization  of  higher  order  modes  for 
fundamental  mode  discrimination. 


7.  Conclusion 

A  model  of  transient  modal  instability  in  fiber  amplifiers  has  been  presented.  This  model 
relies  on  a  time-dependent  3-dimensional  treatment  of  the  effects  of  quantum  defect  heating 
on  the  waveguiding  properties  of  amplifier  fibers.  This  model  has  confirmed  some 
predictions  first  made  using  coupled-mode  theory.  These  include  the  existence  of  stable  and 
unstable  regions  along  the  length  of  the  fiber,  their  evolution  over  time,  the  increase  in  the 
severity  of  the  instability  toward  the  output  end  of  the  fiber*  and  rapid  optical  intensity 
variations  along  the  fiber.  It  has  also  exhibited  several  additional  aspects  of  TMI  including 
coupling  to  un-seeded  modes  and  lower  onset  threshold  powers  compared  to  a  prior  model 
that  are  more  in  line  with  experimental  observations. 

Efforts  have  been  made  to  reduce  the  time  required  to  perform  the  calculations  including 
the  realization  of  a  new  hybrid  finite  element,  harmonic  beam  propagation  method,  a  hybrid 
finite  element,  harmonic,  finite  difference  thermal  solver,  and  parallel  implementation  on 
modern  high  performance  computing  architectures  using  the  message  passing  interface.  This 
method  should  prove  a  useful  tool  in  studying  and  eventually  overcoming  the  technological 
challenges  presented  by  TMI  in  high  average  power  fiber  amplifiers. 
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